################################################################################
###
### Filename: 	deeper-roots-ABS-replication.R
### Author: 	  David A. Bateman
### Function: 	Replicates and extends analyses in Political Legacy of American 
###				      Slavery (ABS 2016) for use in online appendix of Deeper Roots   
###             by Bateman and Schickler.
###             Most of this code is based on slavery-jop-replication.R, 
###             submitted by Acharya, Blackwell, and Sen to the Journal of 
###             Politics dataverse to reproduce figures and tables in Acharya, 
###             Blackwell, and Sen (2016). 
###	Last update: January 1, 2023
###
################################################################################

rm(list = ls())

library(ggpubr)
library(foreign)
library(plyr)
library(reshape)
library(sandwich)
library(maps)
library(stargazer)
library(AER)
library(Formula)
library(lme4)
library(cem)
library(latticeExtra)
library(stringr)
data(state.fips)
state.fips <- unique(state.fips[,c("fips","abb")])
state.fips$abb <- as.character(state.fips$abb)
state.fips <- rbind(state.fips, c(2, "AK"))
state.fips <- rbind(state.fips, c(15, "HI"))
rownames(state.fips) <- state.fips$abb
fips.state <- state.fips
rownames(fips.state) <- fips.state$fips
data(county.fips)
dodgerblue.30 <- rgb(30, 144, 255, 76.5, max =255)
indianred.30 <- rgb(205, 92, 92, 76.5, max =255)
indianred.75 <- rgb(205, 92, 92, 191, max =255)

dir <-paste("../Deeper Roots/", sep="")
source("Dataverse/ABS-Files/panel-utils.R")

## http://people.su.se/~ma/clustering.pdf
robust.se <- function(fm, clvar){
  # R-codes (www.r-project.org) for computing
  # clustered-standard errors. Mahmood Arai, Jan 26, 2008.
  # The arguments of the function are:
  # fitted model, cluster1 and cluster2
  # You need to install libraries `sandwich' and `lmtest'
  library(sandwich);library(lmtest);
  x <- eval(fm$call$data, envir = parent.frame())
  if ("polr" %in% class(fm)) {
    require(MASS)
    cluster <- x[rownames(predict(fm, type = "probs")), clvar]
  } else {
    cluster <- x[names(predict(fm)), clvar]
  }
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- dim(vcov(fm))[1]
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL)
}

ch.row <- function(name, yesno) {
  c(name, ifelse(yesno, "$\\checkmark$", ""))
}

countydata <- read.csv("Dataverse/ABS-Files/abs-jop-countydata-ext-1830.csv", stringsAsFactors = FALSE)
wh.counties <- read.csv("Dataverse/ABS-Files/abs-jop-cces-white-countydata-ext-1830.csv", stringsAsFactors = FALSE)
cces.comb <- read.csv("Dataverse/ABS-Files/abs-jop-cces-ind-ext-1830.csv", stringsAsFactors = FALSE)

st.list <- c("AL", "AR", "GA", "FL", "KY", "LA", "MS", "MO", "NC", "SC", "TN", "TX", "VA","WV")
cces.comb$abs.sample <- 1 * (cces.comb$state.abb %in% st.list)
wh.counties$abs.sample <- 1 * (wh.counties$state.abb %in% st.list)
countydata$abs.sample <- 1 * (countydata$state.abb %in% st.list)
wh.counties$tractor.growth <- (wh.counties$tractors40 - wh.counties$tractors30)

cces.comb$inc.cat <- factor(cces.comb$inc.cat, levels = c("<20k", "20-50k", "50-100k", "100-150k", "150k+"))
whites <- cces.comb[which(cces.comb$white == 1),]
blacks <- cces.comb[which(cces.comb$black == 1),]
latinos <- cces.comb[which(cces.comb$latino == 1),]
others <- cces.comb[which(cces.comb$white != 1 & cces.comb$black != 1 & cces.comb$latino != 1),]

dim(cces.comb)
dim(whites)
sum(whites$state.abb %in% st.list)

dim(blacks)
sum(blacks$state.abb %in% st.list)

dim(latinos)
sum(latinos$state.abb %in% st.list)


## Individual-level data
southerners <- subset(cces.comb, abs.sample == 1)
s.whites <- subset(whites, abs.sample == 1)
s.blacks <- subset(blacks, abs.sample == 1)
s.latinos <- subset(latinos, abs.sample == 1)
s.whites$state.abb <- factor(s.whites$state.abb)
s.blacks$state.abb <- factor(s.blacks$state.abb)
s.latinos$state.abb <- factor(s.latinos$state.abb)

## County-level data
south.counties <- subset(wh.counties, abs.sample == 1)
south.counties$state.abb <- factor(south.counties$state.abb)
south.counties <- south.counties[order(as.numeric(south.counties$fips)),]

nrow(s.whites)
sum(countydata$state.abb %in% st.list)
dim(south.counties)

## NES Results
nes.counties <- read.csv("Dataverse/ABS-Files/abs-jop-nes-white-countydata-ext-1830.csv", stringsAsFactors = FALSE)
nes.counties$abs.sample <- 1 * (nes.counties$state.abb %in% st.list)
nes.comb <- read.csv("Dataverse/ABS-Files/abs-jop-nes-ind.csv", stringsAsFactors = FALSE)
nes.comb$abs.sample <- 1 * (nes.comb$state.abb %in% st.list)
nes.whites <- nes.comb[which(nes.comb$white == 1),]
nes.blacks <- nes.comb[which(nes.comb$black == 1),]

## Individual-level analysis
ns.whites <- subset(nes.whites, abs.sample == 1)
ns.blacks <- subset(nes.blacks, abs.sample == 1)
ns.whites$state.abb <- factor(ns.whites$state.abb)
ns.blacks$state.abb <- factor(ns.blacks$state.abb)

dim(nes.whites)
dim(ns.whites)
dim(nes.blacks)
dim(ns.blacks)
sum(nes.counties$state.abb %in% st.list)

### Deeper Roots: Difference between 1860 and 1830 enslavement
south.counties$diff <- south.counties$pslave1860 - south.counties$pslave1830
nes.counties$diff <- nes.counties$pslave1860 - nes.counties$pslave1830

base1860.form <- formula(. ~ pslave1860 + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2)+ rugged  + land.ineq1860 + sfarmprop1860 + log(totpop1860) + log(fvalpac1860) + log(acimp1860) + fbprop1860  + rail1860 + water1860 + state.abb)
base1860.1830.form <- formula(. ~ pslave1860 + pslave1830 + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2)+ rugged  + land.ineq1860 + sfarmprop1860 + log(totpop1860) + log(fvalpac1860) + log(acimp1860) + fbprop1860  + rail1860 + water1860 + state.abb)
base1860.diff.form <- formula(. ~ diff + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2)+ rugged  + land.ineq1860 + sfarmprop1860 + log(totpop1860) + log(fvalpac1860) + log(acimp1860) + fbprop1860  + rail1860 + water1860 + state.abb)

ind.form <- formula(. ~ pslave1860   + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2) + rugged + land.ineq1860 + sfarmprop1860 + log(totpop1860) + log(fvalpac1860) + log(acimp1860) + fbprop1860 + rail1860 + water1860 + as.factor(educ) +  inc.cat +religion + female + age + state.abb*as.factor(year))

ind.int.form <- formula(. ~ pslave1860   + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2) + rugged + land.ineq1860 + sfarmprop1860 + log(totpop1860) + log(fvalpac1860) + log(acimp1860) + fbprop1860 + rail1860 + water1860 + as.factor(educ)*pslave1860 + inc.cat*pslave1860 + religion*pslave1860 + female*pslave1860 + age*pslave1860  + state.abb*as.factor(year))

context.form <- formula(. ~ pslave1860   + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2) + rugged + land.ineq1860 + sfarmprop1860 + log(totpop1860) + log(fvalpac1860) + log(acimp1860) + fbprop1860 + rail1860 + water1860 + as.factor(educ) + inc.cat  +religion +female + age + blkprop.z00 + log(medinc.z10) + w.unemp.rate2014 + log(wbincratio2014) + state.abb*as.factor(year))
context.int.form <- formula(. ~ pslave1860   + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2) + rugged + land.ineq1860 + sfarmprop1860 + log(totpop1860) + log(fvalpac1860) + log(acimp1860) + fbprop1860 + rail1860 + water1860 + as.factor(educ) +  inc.cat  +religion +female + age + blkprop.z00*pslave1860 + log(medinc.z10)*pslave1860 + w.unemp.rate2014*pslave1860 + log(wbincratio2014)*pslave1860 + state.abb*as.factor(year))

## have to use Formula package for ivreg calls
base.iv.form <- Formula(. ~ pslave1860 + log(coarea00) + rugged + latitude + I(latitude^2) + longitude + I(longitude^2)  + water1860  + state.abb | cottonsuit + log(coarea00) + rugged  + latitude + I(latitude^2) + longitude + I(longitude^2) + water1860  + state.abb)

### Deeper Roots: Additional models including 1830, difference between 1860-1830
### and 1860 plus this difference
base.iv.form.1830 <- Formula(. ~ pslave1860 + pslave1830 + log(coarea00) + rugged + latitude + I(latitude^2) + longitude + I(longitude^2)  + water1860  + state.abb | cottonsuit + log(coarea00) + rugged  + latitude + I(latitude^2) + longitude + I(longitude^2) + water1860  + state.abb)
base.iv.form.diff <- Formula(. ~ diff + log(coarea00) + rugged + latitude + I(latitude^2) + longitude + I(longitude^2)  + water1860  + state.abb | cottonsuit + log(coarea00) + rugged  + latitude + I(latitude^2) + longitude + I(longitude^2) + water1860  + state.abb)
base.iv.form.diff2 <- Formula(. ~ pslave1860 + diff + log(coarea00) + rugged + latitude + I(latitude^2) + longitude + I(longitude^2)  + water1860  + state.abb | cottonsuit + log(coarea00) + rugged  + latitude + I(latitude^2) + longitude + I(longitude^2) + water1860  + state.abb)

base.first.form <- formula(pslave1860 ~ cottonsuit + log(coarea00) + rugged + latitude + I(latitude^2) + longitude + I(longitude^2)  +water1860 + state.abb)
rform.form <- formula(. ~  cottonsuit + log(coarea00) + rugged + latitude + I(latitude^2)+ longitude + I(longitude^2)  + water1860+  state.abb)
rform.form.1830 <- formula(. ~  cottonsuit + pslave1830 + log(coarea00) + rugged + latitude + I(latitude^2)+ longitude + I(longitude^2)  + water1860+  state.abb)
rform.form.diff <- formula(. ~  cottonsuit + diff + log(coarea00) + rugged + latitude + I(latitude^2)+ longitude + I(longitude^2)  + water1860+  state.abb)

## Tables
cnty.res <- lm(dem ~ pslave1860, data = south.counties, weights = sample.size)
summary(cnty.res)

cnty.res.fe <- lm(dem ~ pslave1860 + state.abb, data = south.counties, weights = sample.size)
summary(cnty.res.fe)

cnty.res.full <- lm(update(base1860.form, dem ~ .), data = south.counties, weights = sample.size)
summary(cnty.res.full)

cnty.aff <- lm(affirm ~ pslave1860, data = south.counties, weights = sample.size)
summary(cnty.aff)

cnty.aff.fe <- lm(affirm ~ pslave1860 + state.abb, data = south.counties, weights = sample.size)
summary(cnty.aff.fe)

cnty.aff.full <- lm(update(base1860.form, affirm ~ .), data = south.counties, weights = sample.size)
summary(cnty.aff.full)

cnty.resent <- lm(resent ~ pslave1860, data = south.counties, weights = sample.size.res)
summary(cnty.resent)

cnty.resent.fe <- lm(resent ~ pslave1860 +  state.abb, data = south.counties, weights = sample.size.res)
summary(cnty.resent.fe)

cnty.resent.full <- lm(update(base1860.form, resent ~ .), data = south.counties, weights = sample.size.res)
summary(cnty.resent.full)

## NES Individual Results
therm.mod <- lm(therm.diff ~ pslave1860, data = ns.whites, weights = weight)
therm.mod.rse <- robust.se(therm.mod, clvar = "fips")
therm.mod.rse

therm.mod.fe <- lm(therm.diff ~ pslave1860 + state.abb*as.factor(year), data = ns.whites, weights = weight)
therm.mod.fe.rse <- robust.se(therm.mod.fe, clvar = "fips")
therm.mod.fe.rse

therm.1860 <- lm(update(base1860.form, therm.diff ~ . + state.abb*as.factor(year)), data = ns.whites, weights = weight)
therm.1860.rse <- robust.se(therm.1860, clvar = "fips")
therm.1860.rse


## for comparison with coefficient
0.2*abs(coef(cnty.res.full)["pslave1860"])/sd(south.counties$dem, na.rm = TRUE)
0.2*abs(coef(cnty.aff.full)["pslave1860"])/sd(south.counties$affirm, na.rm = TRUE)
0.2*coef(cnty.resent.full)["pslave1860"]/sd(south.counties$resent, na.rm = TRUE)
0.2*coef(therm.1860)["pslave1860"]/sd(ns.whites$therm.diff, na.rm = TRUE)

## Table 1 from The Political Legacy of American Slavery (2016)
tab1 <- stargazer(cnty.res, cnty.res.full, cnty.aff.full, cnty.resent.full,
                  keep = "pslave1860", style = "apsr", omit.stat = c("adj.rsq","ll", "F", "ser"),
                  covariate.labels = c("Prop. Slave, 1860"), dep.var.labels = c("Prop Democrat", "Affirm. Action", "Racial Resentment"), column.sep.width = "0pt", float = FALSE, header = FALSE, add.lines = list(rep("", 4), ch.row("State Fixed Effects", c(FALSE, rep(TRUE,3))), ch.row("1860 Covariates", c(FALSE, rep(TRUE,3))), rep("", 4)), multicolumn = TRUE)
tab1[4] <- "\\\\[-1.8ex] & \\multicolumn{2}{c}{Prop. Democrat} & Affirm. Action & Racial Resentment \\\\ "
tab1 <- gsub("\\{\\*\\}", "\\{\\\\dagger\\}", tab1)
tab1 <- gsub("\\{\\*\\*\\}", "\\{\\*\\}", tab1)
tab1 <- gsub("\\{\\*\\*\\*\\}", "\\{\\*\\*\\}", tab1)
tab1 <- tab1[-(length(tab1)-1)]
cat(paste(tab1, collapse = "\n"), "\n")


## Instrumental Variables
## Democratic voting - form in Acharya, Blackwell, Sen (2016) 
cnty.iv <- ivreg(update(base.iv.form, dem ~ .), data = south.counties, weights = sample.size)
summary(cnty.iv)
## Democratic voting - online appendix in Deeper Roots (Bateman and Schickler) 
cnty.iv.1830 <- ivreg(update(base.iv.form.1830, dem ~ .), data = south.counties, weights = sample.size)
summary(cnty.iv.1830)
cnty.iv.diff <- ivreg(update(base.iv.form.diff, dem ~ .), data = south.counties, weights = sample.size)
summary(cnty.iv.diff)
cnty.iv.test <- ivreg(update(base.iv.form.diff2, dem ~ .), data = south.counties, weights = sample.size)
summary(cnty.iv.test)

## Affirmative Action - form in Acharya, Blackwell, Sen (2016)
cnty.aff.iv <- ivreg(update(base.iv.form, affirm  ~ .), data = south.counties, weights = sample.size)
summary(cnty.aff.iv)
## Affirmative Action - online appendix in Deeper Roots (Bateman and Schickler) 
cnty.aff.1830.iv <- ivreg(update(base.iv.form.1830, affirm  ~ .), data = south.counties, weights = sample.size)
summary(cnty.aff.1830.iv)
cnty.aff.diff.iv <- ivreg(update(base.iv.form.diff, affirm  ~ .), data = south.counties, weights = sample.size)
summary(cnty.aff.diff.iv)
cnty.aff.test.iv <- ivreg(update(base.iv.form.diff2, affirm  ~ .), data = south.counties, weights = sample.size)
summary(cnty.aff.test.iv)

## Racial Resentment - form in Acharya, Blackwell, Sen (2016) 
cnty.resent.iv <- ivreg(update(base.iv.form, resent  ~ .), data = south.counties, weights = sample.size.res)
summary(cnty.resent.iv)
## Racial Resentment - online appendix in Deeper Roots (Bateman and Schickler)
cnty.resent.1830.iv <- ivreg(update(base.iv.form.1830, resent  ~ .), data = south.counties, weights = sample.size.res)
summary(cnty.resent.1830.iv)
cnty.resent.diff.iv <- ivreg(update(base.iv.form.diff, resent  ~ .), data = south.counties, weights = sample.size.res)
summary(cnty.resent.diff.iv)
cnty.resent.test.iv <- ivreg(update(base.iv.form.diff2, resent  ~ .), data = south.counties, weights = sample.size.res)
summary(cnty.resent.test.iv)

## Feeling thermometer - form in Acharya, Blackwell, Sen (2016)
therm.cty.iv <- ivreg(update(base.iv.form, therm.diff ~ .), data = nes.counties, weights = sample.size.td, subset = abs.sample == 1)
summary(therm.cty.iv)
## Feeling thermometer - online appendix in Deeper Roots (Bateman and Schickler)
therm.cty.1830.iv <- ivreg(update(base.iv.form.1830, therm.diff ~ .), data = nes.counties, weights = sample.size.td, subset = abs.sample == 1)
summary(therm.cty.1830.iv)
therm.cty.diff.iv <- ivreg(update(base.iv.form.diff, therm.diff ~ .), data = nes.counties, weights = sample.size.td, subset = abs.sample == 1)
summary(therm.cty.diff.iv)
therm.cty.test.iv <- ivreg(update(base.iv.form.diff2, therm.diff ~ .), data = nes.counties, weights = sample.size.td, subset = abs.sample == 1)
summary(therm.cty.test.iv)


cnty.res.first <- lm(base.first.form, data = south.counties)
summary(cnty.res.first)

################################################################################
### Figure A18: Acharya, Blackwell, and Sen’s instrumented coefficient       ###
### estimates for 1860 slavery (2016, Table 2)                               ###
################################################################################

x<-coefplot::multiplot(cnty.iv,cnty.aff.iv,cnty.resent.iv, by = c("Model"),lwdInner = 1, lwdOuter = 1, intercept=FALSE, ylab="", coefficients = c("pslave1860"), title = "", names=c("Prop. Dem","Aff. Action", "Racial Resent"))
y<-coefplot::coefplot(therm.cty.iv, by = c("Model"), lwdInner = 1, lwdOuter = 1,intercept=FALSE, coefficients = c("pslave1860"), ylab="", title = "", newNames=c(pslave1860="Diff. Thermo."))
plot1<- ggarrange(x,y, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
cairo_pdf(file = "Dataverse/Figures/Figure-A18.pdf", width = 8, height = 4, pointsize = 8, family = "Minion Pro")
annotate_figure(plot1, top = text_grob("Acharya, Blackwell, and Sen's instrumented coefficient estimates for 1860 slavery (2016, Table 2)",size = 10))
ggarrange(x,y, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
dev.off()

################################################################################
### Figure A19: Instrumented estimates when controlling for earlier          ###
### patterns of enslavement                                                  ###
################################################################################

a<-coefplot::multiplot(cnty.iv.1830,cnty.aff.1830.iv,cnty.resent.1830.iv,outerCI=1.96,lwdInner = 1, lwdOuter = 1, by = c("Model"), intercept=FALSE, ylab="Row (1)", coefficients = c("pslave1860"), title = "", names=c("Prop. Dem","Aff. Action", "Racial Resent"))
b<-coefplot::coefplot(therm.cty.1830.iv, by = c("Model"),outerCI=1.96,lwdInner = 1, lwdOuter = 1, intercept=FALSE, ,coefficients = c("pslave1860"), ylab="", title = "", newNames=c(pslave1860="Diff. Thermo."))
plot2<- ggarrange(a,b, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
annotate_figure(plot2, top = text_grob("Instrumented coefficient estimates for 1860 slavery, controlling for 1830 enslavement",size = 12))

c<-coefplot::multiplot(cnty.iv.diff,cnty.aff.diff.iv,cnty.resent.diff.iv,by = c("Model"), lwdInner = 1, lwdOuter = 1,intercept=FALSE, ylab="Row (2)", coefficients = c("diff"), title = "", names=c("Prop. Dem","Aff. Action", "Racial Resent"))
d<-coefplot::coefplot(therm.cty.diff.iv, by = c("Model"),outerCI=1.96,lwdInner = 1, lwdOuter = 1, intercept=FALSE, coefficients = c("diff"), ylab="", title = "", newNames=c(diff="Diff. Thermo."))
plot3<- ggarrange(c,d, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
annotate_figure(plot3, top = text_grob("Instrumented coefficient estimates for change between 1860 and 1830",size = 12))

e<-coefplot::multiplot(cnty.iv.test,cnty.aff.test.iv,cnty.resent.test.iv,by = c("Model"),lwdInner = 1, lwdOuter = 1, intercept=FALSE, ylab="Row (3)", coefficients = c("pslave1860"), title = "", names=c("Prop. Dem","Aff. Action", "Racial Resent"))
f<-coefplot::coefplot(therm.cty.test.iv, by = c("Model"),outerCI=1.96,lwdInner = 1, lwdOuter = 1, intercept=FALSE, coefficients = c("pslave1860"), ylab="", title = "", newNames=c(pslave1860="Diff. Thermo."))
plot4<- ggarrange(e,f, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
annotate_figure(plot4, top = text_grob("Instrumented coefficient estimates for 1860 slavery, controlling for change from 1830",size = 12))


cairo_pdf(file = "Dataverse/Figures/Figure-A19.pdf", width = 8, height = 11, pointsize = 8, family = "Minion Pro")
ggarrange(plot2, plot3, plot4, ncol=1)
dev.off()

################################################################################
## The following reproduces ABS (2016) Figure 4, for the full South and different subsets. 
## The code is from slavery-jop-replication.R and modified to subset different states. 

## 1 = all southern states
## 2 = just those states examined in "Deeper Roots"
## 3 = just those states for which we have 1850s voting data
## 4 = just those states for which we have roll call votes on Free Black rights / slavery

for (i in 1:6) {
  vers = i
  if (vers==1) {
    st.list1 <- c("AL", "AR", "GA", "FL", "KY", "LA", "MS", "MO", "NC", "SC", "TN", "TX", "VA","WV")
    title <- "All Southern States"
    file_name <- "Dataverse/Figures/Figure-A20.pdf"
    year.list <- seq(1840, 1964, by = 4)
  }
  if (vers==2) {
    st.list1 <- c("AL", "GA", "KY", "MS", "NC", "SC", "TN", "VA","WV")
    title <- "All States Analyzed in 'Deeper Roots'"
    file_name <- "Dataverse/Figures/Figure-A21.pdf"
    year.list <- seq(1840, 1964, by = 4)
  }
  if (vers==3) {
    st.list1 <- c("AL", "GA", "KY", "MS", "TN", "SC")
    title <- "Southern Rights Elections"
    file_name <- "Dataverse/Figures/Figure-A23.pdf"
    year.list <- seq(1840, 1964, by = 4)
  }
  if (vers==4) {
    st.list1 <- c("NC", "TN", "VA","WV", "MS")
    title <- "Roll Call Vote States"
    file_name <- "Dataverse/Figures/Figure-A22.pdf"
    year.list <- seq(1840, 1964, by = 4)
  }

  ## Most county boundaries are set by 1925 or so.
  ## variables without the 1860 suffix are based on the 1860 county
  ## boundaries

  pres.form <- formula(. ~ sprop + log(coarea00) + latitude + I(latitude^2) + longitude + I(longitude^2)+ rugged  + land.ineq + sfarmprop + log(totpop) + log(fvalpc) + log(acimp) + fbprop  + rail + water + state.abb)
  pres.iv.form <- Formula(. ~ sprop + log(coarea00) + rugged + latitude + I(latitude^2) + longitude + I(longitude^2)  + water  + state.abb | cottonsuit + log(coarea00) + rugged  + latitude + I(latitude^2) + longitude + I(longitude^2) + water  + state.abb)
  outvars <- paste("pdem", year.list, sep = "")
  pdemcoefs <- matrix(NA, nrow = length(outvars), ncol = 3)
  pdemcoefs.nox <- matrix(NA, nrow = length(outvars), ncol = 3)
  pdemcoefs.iv <- matrix(NA, nrow = length(outvars), ncol = 3)
  pdemcoefs.rfns <- matrix(NA, nrow = length(outvars), ncol = 3)
  pdemcoefs.rf <- matrix(NA, nrow = length(outvars), ncol = 3)
  for (y in 1:length(outvars)) {
    if (!(outvars[y] %in% c("pdem1864", "pdem1868"))) {
      ## OLS
      ff <- as.formula(paste(outvars[y], " ~ ."))
      if (year.list[y] < 1924) {
        thismod <- lm(update(pres.form, ff), data = countydata, subset = state.abb %in% st.list1)
      } else {
        thismod <- lm(update(base1860.form, ff), data = countydata, subset = state.abb %in% st.list1)
      }
      pdemcoefs[y,1] <- 0.25*coef(thismod)[2]
      pdemcoefs[y,2:3] <- 0.25*confint(thismod)[2,]
      ## Only state FEs
      if (year.list[y] < 1924) {
        ff <- as.formula(paste(outvars[y], " ~ sprop + state.abb"))
      } else {
        ff <- as.formula(paste(outvars[y], " ~ pslave1860 + state.abb"))
      }
      thismod <- lm(ff, data = countydata, subset = state.abb %in% st.list1)
      pdemcoefs.nox[y,1] <- 0.25*coef(thismod)[2]
      pdemcoefs.nox[y,2:3] <- 0.25*confint(thismod)[2,]
      ## IV
      if (year.list[y] < 1924) {
        ff <- update(pres.iv.form, as.formula(paste(outvars[y], "~ . | .")))
      } else {
        ff <- update(base.iv.form, as.formula(paste(outvars[y], "~ . | .")))
      }
      thismod <- ivreg(ff, data = countydata, subset = state.abb %in% st.list1)
      pdemcoefs.iv[y,1] <- coef(thismod)[2]
      pdemcoefs.iv[y,2:3] <- confint(thismod)[2,]
      ## Reduced form in the non-defined part of the South
      ff <- as.formula(paste(outvars[y], "~ cottonsuit + log(coarea00) + rugged+ latitude + I(latitude^2) + longitude + I(longitude^2)   + state.abb"))
      if (vers==1) {
        thismod <- lm(ff, data = countydata, subset = !(state.abb %in% st.list1))
      }
      if (vers!=1) {
        thismod <- lm(ff, data = countydata, subset = !(state.abb %in% st.list1) & (state.abb %in% st.list))      
      }
      pdemcoefs.rfns[y,1] <- coef(thismod)[2]
      pdemcoefs.rfns[y,2:3] <- confint(thismod)[2,]
      ## Reduced form in the defined part of the South
      thismod <- lm(ff, data = countydata, subset = state.abb %in% st.list1)
      pdemcoefs.rf[y,1] <- coef(thismod)[2]
      pdemcoefs.rf[y,2:3] <- confint(thismod)[2,]
    }
    ## if (outvars[y] %in% c("pdem1860","pdem1864")) {
    ##   pdemcoefs[y,1] <- NA
    ##   pdemcoefs[y,2:3] <- NA
    ##   pdemcoefs.iv[y,1] <- NA
    ##   pdemcoefs.iv[y,2:3] <- NA
    ## }
    ## if (outvars[y] == "pdem1856") {
    ##   pdemcoefs[y,1] <- NA
    ##   pdemcoefs[y,2:3] <- NA
    ## }
  }
  
  douglas.ols <- lm(update(pres.form, pdem1860 ~ .), data = countydata, subset = state.abb %in% st.list1)
  douglas.iv <- ivreg(update(pres.iv.form, pdem1860 ~ .), data = countydata, subset = state.abb %in% st.list1)
  
  wallace.ols <- lm(update(base1860.form, wallace68.alt ~ .), data = countydata, subset = state.abb %in% st.list1)
  wallace.iv <- ivreg(update(base.iv.form, wallace68.alt ~ .), data = countydata, subset = state.abb %in% st.list1)
  thurmond.ols <- lm(update(base1860.form, thurmond48 ~ .), data = countydata, subset = state.abb %in% st.list1)
  thurmond.iv <- ivreg(update(base.iv.form, thurmond48 ~ .), data = countydata, subset = state.abb %in% st.list1)
  obama.ols <- lm(update(base1860.form,  wht.obama.vote ~ .), data = countydata, subset = abs.sample == 1)
  obama.iv <- ivreg(update(base.iv.form, wht.obama.vote  ~ .), data = countydata, subset = abs.sample == 1)
  
  cairo_pdf(file= file_name, family = "Minion Pro", height = 4.5, width = 6.5, pointsize = 9)
  plot(x = year.list, y = 0.25*pdemcoefs.iv[,1], ylim = range(c(.25*pdemcoefs.iv,25*confint(obama.iv)["pslave1860",]), na.rm = TRUE), xlim=c(min(year.list),2016), xlab = "Year",
       ylab = "Effect of Slavery on % Democrat", pch = 19, main = title, bty = "n", yaxt = "n")
  abline(v = 1904, lty = 2, col = "grey70")
  abline(v = 1965, lty = 2, col = "grey70")
  axis(side = 2, las = 2, cex = 0.8)
  abline(h=0, col = "grey")
  segments(x0 = year.list, y0 = .25*pdemcoefs.iv[,2], y1 = .25*pdemcoefs.iv[,3])
  rect(xleft = 1860, xright = 1877, ybottom = -100, ytop=100, col = rgb(.5,.5,.5, alpha = 0.5), border = NA)
  points(x = 1968, y = .25*coef(wallace.iv)["pslave1860"], pch = 17, col = "indianred")
  segments(x0 = 1968, y0 = 0.25*confint(wallace.iv)["pslave1860",1], y1 = 0.25*confint(wallace.iv)["pslave1860",2], col = "indianred")
  text(x = 1968, y = 0.25*coef(wallace.iv)["pslave1860"], "Wallace\n1968", pos = 4, col = "indianred")
  points(x = 1949, y = .25*coef(thurmond.iv)["pslave1860"], pch = 17, col = "indianred")
  segments(x0 = 1949, y0 = 0.25*confint(thurmond.iv)["pslave1860",1], y1 = 0.25*confint(thurmond.iv)["pslave1860",2], col = "indianred")
  text(x = 1949, y = 0.35*confint(thurmond.iv)["pslave1860",2], "Thurmond\n1948", pos = 3, col = "indianred")
  segments(x0=1949, y0=0.35*confint(thurmond.iv)["pslave1860",2],y1=0.26*confint(thurmond.iv)["pslave1860",2], lty = 3, col = "grey70")
  points(x = 2008, y = 25*coef(obama.iv)["pslave1860"], pch = 19)
  segments(x0 = 2008, y0 = 25*confint(obama.iv)["pslave1860",1], y1 = 25*confint(obama.iv)["pslave1860",2])
  text(x = 2008, y = 25*coef(obama.iv)["pslave1860"], "Obama\n2008", pos = 2)
  dev.off()
  

}
